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ABSTRACT 


In  designing  vehicles  capable  of  weathering  rough  seas  produced 
by  intense,  persevering  storms,  some  method  of  generating  a repre- 
sentation  of  the  important  effects  of  a random  sea  must  be  employed . 
The  proposed  technique  models  the  wave  amplitude  as  a function  of 
horizontal  location  and  time  and  evaluates  the  accompanying  hori- 
zontal and  vertical  components  of  orbital  particle  velocity  as  a 
function  of  horizontal  and  vertical  location  and  time.  The  under- 
lying definition  of  a random  sea  is  in  terms  of  a Gaussian,  station- 
ary random  process  with  a Pierson-Moskowitz  spectral  density  func- 
tion for  the  wave  amplitude  at  some  reference  point.  The  method  is 
suitable  for  computing  probability  distributions  for  the  vehicle 
state  variables  whether  or  not  the  vehicle  dynamics  are  governed  by 
a linear  differential  equation.  If  the  dynamics  are  linear,  the 
actual  Gaussian  stationary  probability  distribution  can  be  computed 
through  solution  of  an  algebraic  system.  For  the  nonlinear  case, 
random  number  generation  and  numerical  integration  produce  time 
averages  to  approximate  ensemble  averages  in  estimating  any  number 
of  statistical  moments. 


GENERATION  OF  WAVE  AMPLITUDE  AND  ORBITAL 
PARTICLE  VELOCITY  FIELD  OF  A RANDOM  SEA 

Over  the  past  several  years  in  activities  related  to  hydrofoil 
design,  Grumman  engineers  have  amassed  a wealth  of  experience  in  pre- 
dicting the  statistics  of  the  response  of  vehicles  to  random  seas. 

Most  of  this  work  has  been  in  designing  ships  that  use  fully  submerged 
hydrofoils  for  support  (Ref.  1).  Currently,  this  work  is  being  ap- 
plied to  the  design  of  a buoy.  We  have  developed  a unique  mathemati- 
cal representation  of  randan  seas  that  allows  precise  modeling  of  dis- 
tributed effects  on  water-supported  vehicles.  This  representation  is 
far  superior  to  the  usual  filtered  white  noise  models  that  must  rely 
on  such  artifices  as  average  wave  celerity  to  estimate  phase  shift  as 
a function  of  location. 

Our  representation  is  an  adaptation  of  Papoulis* s approximate 
Fourier  expansion  (Ref.  2),  which  leads  to  a mean-square  almost- 
periodic  random  process.  We  obtain  a better  wave  representation  by 
taking  N equiprobable  frequencies  instead  of  N equally  spaced 
ones.  In  fact,  Gikhman  and  Skorokhod  (Ref.  3)  show  that  any  such 
scheme  of  choosing  frequencies  leads  to  a convergent  approximation 
to  the  random  process  as  N increases  without  bound;  but  in  prac- 
tice, the  equiprobable  method  converges  faster  as  a function  of  N. 
Figure  1 represents  a typical  convergence  rate  of  the  truncation 
error  in  vehicle  response  statistics. 

Perhaps  the  best  way  to  describe  the  determination  of  the  ap- 
proximate random  process  is  to  develop  graphically  a fifth  order 
representation.  The  actual  19-knot  Pierson-Moskowitz  spectral 
density  (Ref.  4)  representing  a typical  sea  state  4 is  depicted 
in  Fig.  2a.  Figure  2b  is  the  corresponding  spectral  function  (the 
integral  of  the  spectral  density  in  the  notation  of  Gikhman  and 
Skorokhod).  The  ordinate  of  the  spectral  function  is  divided  into 
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Fig.  1 Relative  Error  in  the  Response  Statistics 
versus  Order  of  the  Approximation 

five  equal  parts  from  its  origin  to  its  final  value  (which  is  the 
2 

variance,  o , of  the  process).  The  center  of  each  of  these  par- 
titions along  the  ordinate  determines  five  discrete  frequencies, 
ox.  (The  dotted  lines  on  the  figure  shows  how  these  are  determined.) 
The  resulting  approximate  spectral  function  is  a staircase  with  equal 
height  and  variable  width  steps.  To  the  left  of  each  step  the  given 
spectral  function  is  underestimated,  while  to  the  right  it  is  over- 
estimated by  the.  same  amount.  The  approximate  spectral  density  is 
a train  of  equal  strength  impulses  at  o^,  and  op.. 

Note  that  this  method  concentrates  the  frequencies  near  the  peak 
of  the  spectral  density.  The  actual  determination  of  the  cd^'s  for 
any  N has  been  computerized  using  the  Newton-Raphson  algorithm. 

2 

Once  the  variance,  o , the  order,  N,  and  the  frequencies, 
i = 1,  2,  . . . , N,  are  specified,  the  approximate  random  process 
for  wave  amplitude,  r) , as  a function  of  time,  t,  takes  the  form: 
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N 

tj  (t)  = } {a.  cos  o>.  t + b.  sin  CD.t] 

W X XX  X 

i=l 

where  the  a^’s  and  bJs  are  identically  distributed  and  uncor- 
related Gaussian  random  variables  with  zero  mean  and  variance, 

2 

a /N.  This  random  process  is  stationary  with  zero  mean  and  auto- 
correlation: 


N 

R(t)  = (cr^/N)  Y,  cos  (^T) 
i=l 

It  is  also  ergodic  in  the  mean,  i.e.,  the  time  average 
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tj  ( t)  dt 
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a.  sin  cd.T 
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tends  to  the  ensemble  mean,  zero,  as  T approaches  infinity  for 

all  sample  functions,  i.e.,  independent  of  the  random  selection  of 

the  a^'s  and  b. 's.  S has  a Gaussian  distribution  with  mean 

1 2 

zero  and  variance,  Og. 
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sin  cd.T 
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cd^T 
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2 . [2al2 
s < L'Mtt] 

for  T > Mtt/2co  . , M = 1,  3,  5,  ...  . In  particular  for  T = 

“^min’  4 < k/T<W2-  For  N 4 5>  “min  = °’75>  and  ^ = 

2.80.  Taking  M = 5,  T = 10.5,  and  Og  < 0.04515  so  that 
Pr{ | Sj  > 0.1}  < 0.027. 

Furthermore,  in  the  limit  as  N gets  large,  tih is  process  is 
ergodic  in  variance.  This  can  be  shewn  by  forming  the  time  average 
of  the  square  of  rj(t): 


W = ~ T!2(t)dt 
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j r r 2 2 2 2 

( ) a.  cos  oo.t+b.  sin  co.t  + 2a.b.  cos  <jo.  t sin  oo.t 

4,  I l li  l li  i lj 
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+ a.b.  cos  co.t  sin  co.t  > dt 
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If  W = cr  , the  process  would  be  ergodic  in  the  variance.  Note 

co  A 2 

that  based  on  the  distribution  of  a. , b. , W = 2N/a  W has  the 

XI  00 

chi-square  distribution  with  2N  degrees  of  freedom;  therefore, 

2 4 

the  mean  of  W is  a and  its  variance  is  a /N.  This  proves 

oo  „ 

that  lim  W = a . Furthermore,  the  chi-square  distribution  sup- 
N-+eo  “ 

plies  a quantitative  estimate  of  the  necessary  order,  N,  to 
guarantee  "practical"  ergodicity  in  variance,  given  that  the  a^'s 
and  b^'s  are  chosen  completely  at  random.  Suppose  that  a 25  per- 
cent error  can  be  tolerated  if  it  occurs  no  more  than  35  percent 
of  v ,e  time:  0.65  ^ Pr(0.75  a2  < W < 1.25  a2}  = Pr(1.5  N < W < 2.5  N] . 

CO 

Thus  N = 15  would  suffice.  Also,  if  N is  taken  to  be  50,  then 
the  same  25  percent  error  would  only  occur  8 percent  of  the  time, 
and  a 10  percent  error  would  occur  with  a 50/50  probability. 

Of  course  for  computer  simulation  of  random  seas,  a pseudo-random 
sampling  algorithm  can  "stack  the  deck"  by  rejecting  unfavorable 
sets  of  (a^,b^) ' s . 
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To  apply  the  random  process  to  simulation  of  buoy  motion, 
it  is  best  to  put  the  process  -q(t)  into  state  variable  form. 

In  the  following,  the  great  power  of  this  method  will  be  seen  in 
its  ability  to  represent  exactly  the  distributed  effects  of  a 
random  sea.  Water  particle  velocity  is  needed  at  each  of  13 
stations  along  the  length  of  the  buoy.  In  addition,  the  buoy 
experiences  pitch,  heave,  and  surge.  What  is  needed  are  expres- 
sions for  wave  amplitude,  t]  (t),  at  an  arbitrary  horizontal  iner- 
tial  distance,  x,  from  the  reference  and  for  vertical  and  hori- 
zontal components  of  orbital  particle  velocity,  W^  (t)  and 

U-  (t) , respectively,  at  an  arbitrary  horizontal  distance,  x, 

*xz 

from  the  reference  and  vertical  distance,  z,  below  the  mean  water 
surface. 


Consider  the  auxiliary  linear  system 


V » 9M  ; V(0)  = V 


where 


n « 


Cl, 


Clr 


0 


0 


and  each  Cl.  is 

l 


ni  = 


0 a). 


-to.  0 

i 


If  the  initial  condition  vector,  V°,  has  zero  mean  and  co- 
2 

variance  matrix  (a  /N)I,  then  the  random  process  for  wave  ampli- 
tude at  the  reference  point  is: 


7 


Recall  that  the  relationship  linking  wavelength,  A,  to  frequency, 

2 


0),  is 


* = 2g /a> 

where  g is  the  acceleration  of  gravity.  Therefore  the  wave  ampli 
tude  at  any  x is 

" I V°i-1  cos(“it  + + V2i  sltl(V  + ^ ) 


i=l 


N 


- I «“(¥)  V2i-1«>  + M¥)  v2i(0 

x i 


i=l 


The  vertical  component  of  orbital  velocity  at  any  (x,z)  is 

w (t)  - l V2i-1  sinKt+ • V2i  cos(v  + 2f:) 
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Sln(  Aj)  V2i-1^  ’ C0S(  A_.)  V2i^ 


Similarly,  the  horizontal  component  is 
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su 
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S3 


u (t) 

T) 

'xz 


■ 1 exK^K  v2i-i  cosKt  + * v2i  sin(v  + ^) 


= l exp (^K  - c»<^)  v2i-l(t)  + sin(^r)  V2i(t) 


For  the  sake  of  discussion,  suppose  that  20  z-locations  and 
5 xrlocations  were  required.  Then  the  above  three  equations  would 
define  an  output  vector  with  205  components,  i.e., 


W - T1  ,T|  , « • o , T)  ,VJ  ,W  , • • • , W 

'x,  'X0’  'X,  T)  T)  ’ T) 

12  5 'x.,z,  'x,  z_  1 


...,  u 


"1-1  12 


X5Z20  xlzl 


W =HV 

where  H is  a 205  by  2N  matrix  of  constants. 

For  the  case  in  which  the  buoy  is  considered  to  satisfy  linear 
small  perturbation  differential  equations  (10  in  our  current  model) , 
the  situation  is  vastly  simplified.  Only  the  stationarity  of  the 
random  process  needs  to  be  utilized.  The  covariance  matrix  of  buoy 
state  variables  is  time  invariant,  and  these  variables  have  zero 
means  and  obey  a multivariate  Gaussian  probability  distribution. 

Thus,  the  probability  distribution  is  completely  determined  by  the 
steady-state  covariance  matrix,  which  is  computed  as  the  routine 
solution  of  an  algebraic  system  of  equations.  There  is  no  need  for 
either  random  number  generation  or  numerical  integration.  The  com- 
puter code  for  this  computation  exists  and  is  well  tested. 

On  the  other  hand,  this  method  of  representing  the  random  sea 
is  also  well  suited  to  the  nonlinear  case.  Naturally,  there  are 


i 


4 


H m 


complications  in  the  computations.  Care  must  be  exercised  in 
selecting  a random  sample  for  the  set  (a^,  b\,  i = 1,  2,  . ..,  N), 
which  corresponds  to  a particular  observed  sea.  Large  enough  N 
and  T must  be  employed.  The  simulation  technique  is  of  the 
Monte  Carlo  type  in  which  ergodicity  of  the  mean  and  variance  ,of 
both  the  driving  process  and  the  buoy  variables  are  assumed.  Even 
with  these  assumptions,  the  use  of  a finite  simulation  time  inter- 
val implies  a certain  leyel  of  error  in  the  estimation  of  the  mean 
vector  and  the  covariance  matrix,,  The  means  of  the  buoy  variables 
are  not  necessarily  zero,  and  their  multivariate  distribution  is 
not  Gaussian  and,  therefore,  not  describable  by  a finite  number  of 
statistical  moments.  Nevertheless,  useful  conclusions  can  still 
be  drawn  about  the  statistics  of  buoy  response  through  time  averages 
of  the  response  and  the  squared  response  to  any  particular  "sample 
wave"  generated  by  the  above  technique. 

In  conclusion,  it  is  expected  that  use  will  be  made  of  both 
linear  models  of  buoy  dynamics  for  which  precise  statistical  con- 
clusions can  be  drawn,  £nd  nonlinear  models,  which  account  especially 
for  both  the  limited  excess  buoyancy  and  the  inability  of  a cable  to 
offer  compressive  resistance. 
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